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ABSTRACT 

A new theory for determining the mass function of cosmic structures is presented. It re- 
lies on a realistic treatment of collapse dynamics. Gravitational collapse is analyzed in 
the Lagrangian perturbative framework. Lagrangian perturbations provide an approx- 
imation of truncated type, i.e. small-scale structure is filtered out. The collapse time 
is suitably defined as the instant at which orbit crossing takes place. The convergence 
of the Lagrangian series in predicting the collapse time of a homogeneous ellipsoid is 
demonstrated; it is also shown that third-order calculations are necessary in predict- 
ing collapse. Then, the Lagrangian prediction, with a correction for quasi-spherical 
perturbations, can be used to determine the collapse time of a homogeneous ellipsoid 
in a fast and precise way. Furthermore, ellipsoidal collapse can be considered as a 
particular truncation of the Lagrangian series. Gaussian fields with scale- free power 
spectra are then considered. The Lagrangian series for the collapse time is found to 
converge when the collapse time is not large. In this case, ellipsoidal collapse gives 
a fast and accurate approximation of the collapse time; spherical collapse is found 
to poorly reproduce the collapse time, even in a statistical sense. Analytical fits of 
the distribution functions of the inverse collapse times, as predicted by the ellipsoid 
model and by third-order Lagrangian theory, are given. These will be necessary for a 
determination of the mass function, which will be given in paper II. 
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1 INTRODUCTION 

An important outcome of any cosmological model is the 
distribution of the masses of those collapsed clumps which 
form during gravitational collapse; this quantity is usually 
called mass function (hereafter MF), or multiplicity func- 
tion. These structures correspond to the observed galaxies, 
groups and clusters of galaxies. Thanks to recent efforts, 
the masses of these observed cosmic structures can be esti- 
mated, though with large uncertainties and relying on un- 
certain hypothesis. Galaxy cluster masses can be estimated 
by means of three different methods: estimates based on op- 
tical galaxies, such as virial estimates (Biviano et al. 1993; 
see also Bahcall & Cen 1993), X-ray temperatures (Henry 
& Arnaud 1991), and gravitational lensing (see, e.g., Fort & 
Mellier 1994). Every method has its drawbacks: virial esti- 
mates rely on the delicate hypothesis of virial equilibrium; 
X-ray analyses rely on safer gas-dynamical hypotheses but 
are confined to the inner regions of the cluster; lensing es- 
timates directly probe the gravitational potential, but are 
now available only for a few clusters. Group masses are 
much more difficult to estimate, due to the small number 



of galaxies involved and to the uncertain dynamical status 
of groups. Virial mass estimates, corrected for unvirializa- 
tion, are available for a large number of groups (Pisani et 
al. 1992). Galaxy masses are relatively easier to determine, 
due to the more evolved dynamical status and to the large 
number of tracers (stars instead of galaxies) (see Ashman, 
Persic & Salucci 1993). However, mass estimates are mainly 
limited to a few optical radii, while dark-matter halos prob- 
ably extend further. 

The state of the art of the MF theory presents severe 
problems as well. Cosmological structures are the sites of 
highly non- linear dynamics; it is well known that, for generic 
initial conditions, the Newtonian collapse of a general self- 
gravitating system has no known solution in the highly non- 
linear regime. Analytical approximations and N-body simu- 
lations can help to face the problem in an approximate way. 
In this paper I will mainly focus on realistic analytical ap- 
proximations to gravitational collapse, needed to develop a 
MF theory. Another paper (Monaco 1996, hereafter paper 
II) develops the statistical tools necessary to get a MF. 

Despite the intrinsic difficulties, a heuristic solution of 
the MF problem was found as early as 1974 by Press & 
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Schechter (1974; hereafter PS). I have already described in 
full detail the PS MF in a previous paper (Monaco 1995; 
hereafter M95). It suffices now to recall that all the dynam- 
ical difficulties of non-linear dynamics are circumvented by 
assuming that 'something happens' to a mass clump (i.e. 
a fully virialized structure forms) as soon as its linearly 
evolved density contrast S reaches a given threshold S c of or- 
der one. One of the simplest fully non-linear collapse models, 
the spherical model, can be invoked to give a more precise 
determination of the quantity 8 C : indeed, the spherical model 
predicts full collapse to a singularity as soon as the linearly 
extrapolated density contrast reaches 5 C =1.69. Surprisingly, 
the PS formula has been found to reproduce in a more or less 
accurate way the results of large N-body simulations (see ref- 
erences and discussions in M95). Comparisons of the PS MF 
to N-body simulations have often suggested a lower value for 
the S c parameter, even though technical details, such as the 
group-finding algorithm and the shape of the filter function, 
have to be taken into account to obtain a precise number. 
The original PS formulation has been improved in a num- 
ber of papers; in particular, Peacock & Heavens (1990) and 
Bond et al. (1991) have solved the so-called 'cloud-in-cloud' 
problem, while Bower (1991) and Lacey & Cole (1993, 1994) 
have extended the formalism to describing merging histories 
of dark matter halos. A more complete list of references and 
a deep discussion on these points is given in paper II. 

It is worth stressing that the only dynamical ingredi- 
ents of the PS recipe are linear theory and a fixed density 
threshold. Thus, the PS MF can only be considered as a 
fitting formula which nicely describes the results of many 
N-body simulations. A great advantage in using only linear 
theory is the following: as linear theory predicts the density 
field to be rescaled, as it evolves, by a factor which depends 
only on time, the statistical properties of the final field are 
the same as the initial field, which is usually assumed to be a 
Gaussian process. In practice, evolved density fields develop 
strong and non-trivial non-Gaussian features. Nonetheless, 
most works on the MF, since that of PS, have detailed the 
statistical treatment of regions with initial density larger 
than a given threshold (called 'excursion sets'), or peaks of 
the initial density field, rarely detailing the dynamical de- 
scription of collapse. 

Only a few authors have considered more evolved dy- 
namics than linear theory. Lucchin & Matarrese (1988) and 
Porciani et al. (1996), in two different ways, introduced non- 
Gaussianity into the MF. In Lucchin & Matarrese (1988) 
the non-Gaussianity was explicitly introduced in a PS-like 
approach; it could be either primordial or of dynamical ori- 
gin. Porciani et al. (1996) introduced non-Gaussianity in 
the framework of diffusion theory (Bond et al. 1991), by 
putting a reflecting barrier at S — —1, to avoid unphysi- 
cal negative densities. In both cases, gravitationally-induced 
non-Gaussianity was found to cause enhanced production of 
high-mass clumps, and, in Porciani et al. (1996), an intrigu- 
ing cut-off at low masses was predicted. Vergassola et al. 
(1994) calculated the MF in the framework of adhesion the- 
ory; collapsed structures were identified with caustics. An- 
other attempt to model in a more realistic way the formation 
of collapsed clumps was made by Bond & Myers (1996), in 
the framework of peak-patch theory. In that theory, the col- 
lapse times of the peak-patches were estimated by means 
of the homogeneous ellipsoidal collapse model, and the final 



locations were found by means of Zel'dovich (1970) approx- 
imation. Ellipsoidal collapse in a cosmological context was 
used also by Eisenstein & Loeb (1995) to estimate the an- 
gular momentum of collapsing structures. 

A very different approach was used by Cavaliere and 
collaborators in a number of works (see Cavaliere, Menci & 
Tozzi 1994 for a recent review). First they constructed a MF 
theory based on dynamical timescales, then they used a ki- 
netic approach to model the aggregation of already collapsed 
clumps, and finally they found a formalism, based on Cayley 
trees, which was able to unify the kinetic and diffusion ap- 
proaches. They recently applied this formalism to adhesion 
theory (Cavaliere, Menci & Tozzi 1996). Similar or related 
approaches were recently used by Shaviv & Shaviv (1995) 
and by Sheth (1996). Such approaches, very different from 
the ones mentioned above, attempt to describe the highly 
non-linear behaviour of matter clumps after first collapse. 

In a previous paper, M95, I attempted to introduce re- 
alistic dynamics in the framework of a PS-like theory. The 
PS idea of identifying the fraction of collapsed matter as 
the probability that a point fulfills certain collapse condi- 
tions was used to determine the MF. It was shown that 
such kind of approach greatly simplifies if 'local' predictions 
of collapse are used, i.e. if the collapse of a mass element 
is based on some relevant quantities relative only to the el- 
ement considered. Collapse predictions were estimated by 
means of ansatze based on the Zel'dovich approximation and 
by means of the homogeneous ellipsoid collapse model. The 
result was an increase in large-mass objects with respect to 
the usual PS prediction, which could explain why lower <5 C 
values than the spherical 1.69 are found when comparing 
the PS MF to N-body simulations (lower values of S c cor- 
respond to a shift of the large-mass tail of the MF toward 
large masses). 

The M95 work can be seen as a simple 'variation on a 
theme' of the PS procedure, plagued by all the PS faults, 
e.g. the cloud-in-cloud problem. Anyway, that variation suf- 
fices in shedding light on a number of problems which have 
a precise dynamical meaning, and which are connected to 
the use of the spherical collapse model. To be more precise, 
the spherical model predicts full collapse of the whole struc- 
ture to a singularity at a given time; it is usually assumed 
that at that time the structure fully virializes. In this way: 
(i) collapse is quite well defined and (ii) virialization surely 
occurs (iii) at the same moment as the collapse, so that the 
two concepts can be used equivalently. In practice, when col- 
lapse takes place with a realistic geometry, (i) it has to be 
carefully defined, and different authors have in fact differ- 
ent ideas on what collapse is (e.g., is the 'real' collapse of 
a homogeneous ellipsoid that on the first axis, or that on 
all three axes?); (ii) while it is somehow possible to model 
collapse, it is terribly difficult to model virialization, and it 
is terribly difficult to decide when and where virialization 
is going to have place (clusters and groups of galaxies are 
generally not virialized!). 

In this paper, and in paper II, the ideas contained in 
M95 are pushed further, in order to construct a complete 
theory of the mass distributions of cosmic structures, based 
on realistic dynamics. This paper focuses on the dynami- 
cal estimate of the collapse time of a mass element, and on 
the definition of collapse. Paper II will be concerned with 
the statistics necessary to get a MF from the definition of 
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collapse and its distribution. The whole set-up of the dy- 
namical problem is based on the idea, already formulated in 
M95, that the MF is an intrinsically Lagrangian quantity, 
in the sense that it is best faced within the framework of 
Lagrangian fluidodynamics. In Section 2 the machinery of 
Lagrangian dynamics of a cosmological self-gravitating cold 
fluid is presented. Lagrangian perturbation theory is intro- 
duced: it provides an excellent framework to estimate the 
collapse time of a mass element. In Section 3 the defini- 
tion of collapse is analyzed, its punctual, non-local nature 
is stressed, and a way of defining collapse in the case of 
non-filtered fields is presented; this can be implemented as 
the collapse definition when comparing this theory to N- 
body simulations. In Section 4 the Lagrangian perturbative 
series, up to the 3rd order, is successfully applied to the ho- 
mogeneous ellipsoid collapse. In Section 5 the collapse of a 
Gaussian field with power-law spectra, generated in cubic 
grids of 32 3 points, is carefully analyzed. As a result, the 
simple ellipsoidal model is found to give a very accurate es- 
timate of the collapse time, at least for the largest objects 
which collapse. The distribution of the inverses of collapse 
times is analyzed and quantified; this quantity is necessary 
for the calculation of the MF. Section 6 contains a summary 
and conclusions. Technical details about Lagrangian pertur- 
bations and ellipsoidal collapse are given in two Appendices. 



2 THE LAGRANGIAN NATURE OF THE 
MASS FUNCTION 

The MF is the distribution of the masses of the collapsed, 
isolated clumps of mass. To get a MF, information about 
the mass elements which undergo collapse is needed; it is 
not necessary to know where the collapsed clumps go to. 
Then, the MF problem is suitably analyzed in a Lagrangian 
fluidodynamical framework, where the independent space 
variable q is related to a given mass element, whose mo- 
tion is followed. This is different from the usual Eulerian 
picture, where the space variable x is related to a given 
spatial point, and the evolution of the elements, which in- 
stantaneously happen to be in that point, is followed. This 
is the starting point of the MF theory presented here. Re- 
cently, the same point of view has been taken by Audit & 
Alimi (1996). 

Consider a pressureless, unvortical, Newtonian, self- 
gravitating fluid in a perturbed Friedmann-Robertson- 
Walker Universe, with given cosmological parameters flo 
and A. The trajectory x of every (infinitely small) mass el- 
ement, initially at the comoving position q, can be written 
as: 



x(q,t) = q + S(q,t) 



(1) 



S is the displacement field (in comoving coordinates). All 
the kinematic quantities relative to the mass element and 
its density contrast can be expressed i n ter ms of the dis- 
placement S (see Appendix A, equation Al). Note that, in 



this approach, the density is not a dynamical quantity, as it 
was in the Eulerian approach; the only dynamical quantity 
here is the displacement field S. It is possible to find evolu- 
tion equations for S, equivalent to the usual Euler-Poisson 
system for the Newtonian evolution of perturbations. Sev- 
eral authors, namely Buchert (1989), Bouchet et al. (1995), 



Lachieze-Rey (1993b) and Catelan (1995), have written dif- 
ferent equivalent forms of the same system of equations for 
S. The equations of Catelan (1995) are reported in Appendix 
A. In the following it will always be assumed that the initial 
velocity field is irrotational and parallel to the gravitational 
acceleration; see Appendix A for further details. 

It is possible to give initial conditions for the system 
through the initial peculiar rescaled gravitational potential 
y(q, to), instead of the usual initial density contrast <5(q, to); 
the two fields are simply related by a Poisson equation: 

<5(q,*o) 



V ¥>(q,io) = 



b(t ) 



(2) 



where b(t) is the linear growing mode (note that b(to) — 
a(t )). 

An important quantity is the Jacobian determinant of 
the transformation given in equation ([j]): 



■/!q.M =<lrf ( ^ ) =.lot(,\,,. -S.,,,1 



(3) 



(comma denotes q-derivative); the quantity S a ,b is com- 
monly called deformation tensor. The Lagrangian approach 
has a natural limit in the condition J — 0. Before this mo- 
ment, the Lagrangian-to-Eulerian mapping is single-valued 
(single-stream or laminar regime), i.e. mass elements coming 
from different points of (Lagrangian) space do not arrive at 
the same Eulerian position. When J = (orbit crossing or 
shell crossing, hereafter OC) the q — > x mapping becomes 
multi- valued (multi-stream regime), and the density, being 
proportional to the inverse of J, goes to infinity (see Shan- 
darin & Zel'dovich 1989 for a discussion of caustic forma- 
tion). It is beyond the scope of this paper to discuss the pos- 
sibility of using a Lagrangian approach in the multi-stream 
regime; see Buchert (1994) for a discussion. Here it has to 
be noted that OC takes place as soon as the first objects col- 
lapse, which, for any realistic power spectrum, happens very 
soon after recombination. So the Lagrangian formulation is, 
as it stands, essentially useless unless the initial (potential 
or density) field is smoothed to truncate the small-scale part 
of the power spectrum. This is a key point: from a field </3(q), 
which can in principle have power on all scales, a hierarchy 
of smooth fields is worked out: 



v(q) v(q; R 



f) 



(4) 



where Rf is the width of the filter (the two fields are both 
freely called ip, but they are in fact different mathemati- 
cal objects). This highlights the strongest hypothesis of this 
theory, namely that small-scale structure does not influence 
in a significant way the dynamical evolution of larger scales 
before OC. All the following considerations apply to a gen- 
eral element of the ip(q;7?/) hierarchy, i.e. to a smoothed 
version of ip; the parameter Rj will be omitted for simplic- 
ity's sake. (Note that some approximation schemes, such 
as adhesion theory (Gurbatov, Saichev & Shandarin 1989) 
or frozen flow approximation (Matarrese et al. 1992), avoid 
OC; in this case no smoothing of the initial field is needed 
in principle.) 

The Lagrangian evolution equations for S, as well as 
the Eulerian evolution equations for 8, have not been solved 
in general cases. Nonetheless, it is possible to find, as in 
the Eulerian case, perturbative solutions. In the Eulerian 
case (see, e.g., Bouchet 1996), the equations are expanded 
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in terms proportional to powers of the density contrast; if 
this quantity is small, then successive perturbative terms 
are increasingly smaller and the series converges. As a con- 
sequence, the validity of the Eulerian perturbative scheme is 
limited to small density contrasts. In the Lagrangian case, 
the perturbed quantity is not the density, which is not a 
dynamical quantity, but the (comoving) displacements of 
the particles from their initial positions. It is easy to un- 
derstand why this simple change of perturbing parameter 
causes dramatic improvements in the performances of the 
approximations: the de nsity, being proportional to the in- 
verse of J (equation Al), becomes infinite when J goes to 0. 
In this case, at least some matrix element of the deformation 
tensor (e.g. one eigenvalue, if S a ,b is symmetric) is of order 
(minus) one, i.e. just at the limit of validity of the perturba- 
tion scheme. In other words, when Lagrangian perturbations 
start to break down, the density can have reached very large 
values. 

A number of authors have analyzed the Lagrangian 
perturbation scheme at various orders, up to the third 
(Buchert 1989; Moutarde et al. 1991; Bouchet et al 1992; 
Buchert 1992; Buchert & Ehlers 1993; Lachieze-Rey 1993a,b; 
Bouchet et al. 1995; Buchert 1994; Catelan 1995; Bouchet 
1996; Buchert 1996). A brief list of the main results is given 
in Appendix A; see the references for further details. As a 
matter of fact, different authors use very different notations; 
the notation I will use is similar to that of Catelan (1995). 

The perturbative series for S can be written, up to third 
order, as: 

S(q, t) = fax (i)S (1) (q) + b 2 (i)S (2) (q) + h a (t)S (3a) (q) 

+Mi)S (3&) (q) + Mt)S (3c) (q) + • • ■ • (5) 

In the following, among models with A ^ 0, only the flat 
ones will be considered (the others are not of great cosmo- 
logical interest). The time functions b n (t) are given in Ap- 
pendix A, equations (A7). The first-order time function is 
just (minus) the linear growing mode, foi(t) = —b(t), while 
the others are, at leading order and with great accuracy (ex- 
actly for an Einstein-de Sitter background) , proportional to 
b 2 or b 3 , according to their order. This means that, in all 
the calculations that follow, the dependence on the back- 
ground cosmology can be factorized out by using b(t) as 
time variable. The spatial equations for the S^™'(q) terms 
are Poisson equations; this reflects the implicit non-locality 
of Newtonian gravitational dynamics. Again, they are re- 
ported in Appendix A. Note that the S (2) , S (3a) and S (36) 
terms are irrotational is symmetric), while the S*- 3 "' 

term is purely rotational(5'^' is antisymmetric). 

It is easy to recognize that the linear term of the per- 
turbation, 



■b(t)V<p. 



(6) 



is the well-known Zel'dovich (1970) approximation. This is 
not the place to list all the features, merits and limits of 
this approximation; the reviews of Shandarin & Zel'dovich 
(1989) and Sahni & Coles (1995) give full details and com- 
plete reference lists. Some comments were reported also in 
M95. Is is worth mentioning that Zel'dovich and 2nd or- 
der truncated approximations have been found very success- 
ful in predicting the evolution of a N-body matter field in 



the weakly non-linear regime, according to cross-correlation 
tests (Coles, Melott & Shandarin 1993; Melott, Buchert & 
Weifi 1995; Sahni & Coles 1995). The third order has not 
been found to increase significantly the precision of the se- 
ries; moreover, third-order terms are very sensitive to numer- 
ical errors, so the use of third-order predictions is not gen- 
erally recommended in that context. Another general con- 
clusion is that the 3c term does not have great influence on 
the density evolution (see also Buchert et al. 1997), which 
is to be expected, as this term, being purely rotational, cor- 
responds to a rotation of the mass element in Lagrangian 
space, which does not influence the density by itself. 

As already noted in M95, there is a key difference be- 
tween the usual applications of the truncated Lagrangian 
approximations and the one which is to be used here. In the 
works cited above the density fields are evolved up to mass 
variances of order 1 or slightly more; at this level the con- 
vergence of the Lagrangian series is more or less guaranteed 
by construction. But the MF needs a prediction of collapse, 
when all the perturbative terms become of the same order. 
The convergence of the series in this case is not guaranteed 
and has to be checked. This will be done in Section 4 in the 
case of ellipsoidal collapse, and in Section 5 in the case of 
Gaussian fields with scale-free power spectra. 



3 THE DEFINITION OF COLLAPSE 

The key quantity for the determination of a MF is the in- 
stant at which a given mass element collapses. This raises 
the not trivial question of what collapse means. Within the 
Lagrangian perturbation framework, there is a natural def- 
inition of collapse: 



J(q,fe c ) =0 



(7) 



(note that the linear growing mode b is the time variable; b c 
is the collapse time). This instant has already been defined 
before as OC. At this instant a number of things happen: (i) 
the density becomes infinite; (ii) the other kinematic quanti- 
ties become infinite (see M95); (iii) different trajectories in- 
tersect; (iv) multi-stream regions form; (v) shock waves can 
form in a (subdominant) dissipative component (baryons); 
(vi) gravitational dynamics becomes interesting and (vii) re- 
ally difficult to follow. After OC a number of other things can 
happen: violent relaxation, virialization, gas cooling, star 
formation, supernovae feedback etc. All these relevant events 
are decisive in making real astrophysical structures, but are 
very difficult to model. Surely OC is a necessary condition 
for these events to take place. 

The main interest of the present theory is to model 
dark matter clumps, not astrophysical objects. Dark mat- 
ter clumps can be thought of as local high-density concen- 
trations of matter, regardless of their geometry, internal dy- 
namical status (relaxation, virialization or otherwise) and so 
on. Then the above definition of collapse is meaningful, pro- 
vided a mass element which has undergone collapse remains 
in some high density clump; I assume that this is the case. 
This assumption is reasonable: were it not the case, mass 
elements entering a structure would evaporate back into the 
background soon after; such evaporation events are assumed 
to be rare. 
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It is opportune not to introduce conditions on the in- 
ternal dynamical status of a collapsed clump at the level 
of collapse definition, for at least two good reasons: (i) dy- 
namics in the multi-stream regime is not well understood, 
and any oversimplification, such as virialization soon after 
collapse, can be misleading, (ii) Real objects, as groups and 
clusters of galaxies, have plausibly undergone some kind of 
(violent) relaxation, but their state is very complicated; only 
the cores of rich groups and clusters are observed in a more 
evolved state. Furthermore, a collapsed but not dynamically 
evolved clump may contain smaller, more relaxed or fully 
virialized clumps, which can survive in the larger structure 
for a significant time; this could be the situation of galaxies 
in groups. This kind of MF theory cannot provide a descrip- 
tion of such subclumps, or, in other words, cannot represent 
galaxies in groups or clusters. A kinetic approach like that 
described in Cavaliere et al. (1994) can model this kind of 
objects. 

Another key feature of this definition is that it is clearly 
punctual, in the sense that any prediction is relative to a 
point. It is often assumed in the literature that the collapse 
prediction of a point is to be extended to a surrounding re- 
gion of filter- width size, as the prediction is based on a mean 
over a given region; this fact, reasonable especially when 
top-hat smoothing is used, leads to some complications in 
the MF, as Blanchard, Valls-Gabaud & Mamon (1992) and 
Yano, Nagashima & Gouda (1996) have shown. This is not 
my point of view: filtering of the initial field is only necessary 
to suppress small-scale orbit-crossed regions. All the kine- 
matic and dynamical quantities, and then the collapse pre- 
dictions, are strictly punctual, defined on vanishingly small 
mass elements. On the other hand, if the field is smooth on 
a scale Rf, the collapse predictions will show a coherence on 
a scale Rf, without any further assumptions. 

Summing up, the definition of collapse given by equa- 
tion (Q) can reasonably reproduce high-density clumps, re- 
gardless of their internal dynamical status. Due to the trun- 
cated nature of the approximations used, collapsed regions 
are predicted to reach infinite density, but actual collaps- 
ing regions contain smaller-scale structures which spread the 
collapsing mass around, thus lowering the true density. It is 
worthwhile wondering how a clump, which is predicted to 
collapse when the field is smoothed at a given scale, can 
be recognized in the unsmoothed evolved system, e.g., in 
a N-body simulation. Any conventional clump-finding algo- 
rithm, based on percolation or overdensities above a given 
threshold, will probably find more easily those collapsed re- 
gions which have certain geometries, e.g. spherical rather 
than filamentary. The choice of the clump-finding algorithm 
is decisive in this case, and it is opportune to choose an algo- 
rithm based on the collapse definition given above, equation 
([?]). A good choice is to construct an algorithm based sim- 
ply on the concept of mapping, as infinite density is only a 
consequence. I propose the following formula as an imple- 
mentation of equation (^) : 



3qi,q2 : |qi - q»| > L, 

|x(qi,6) - x(q 2 ,6)| < e, 



(8) 



e < L , 



Obviously, a tight correspondence between high-density 
clumps and regions found with equation (^) is expected, es- 
pecially for the most massive clumps. Nonetheless, the pre- 
cise relationship has to be checked with N-body simulations. 
In any case, the 'objects' predicted by this kind of theory are 
interesting in themselves, as they are the sites of highly non- 
linear dynamics, and surely contain the mass which is go- 
ing to form astrophysical and cosmological objects. If more 
restricted classes of objects are needed, the possibility of 
getting them by placing further constraints on the collaps- 
ing point would be worth exploring. As a conclusion, to be 
conservative the MF which can be found with such a defini- 
tion of collapse could be called the MF of high non-linearity 
environments. 



4 LAGRANGIAN PERTURBATIONS IN 

HOMOGENEOUS ELLIPSOIDAL COLLAPSE 

As a first step, the gravitational collapse of a homogeneous 
triaxial ellipsoid is considered. This exercise is useful for test- 
ing the convergence of the Lagrangian series at OC, which 
is not guarenteed by construction. Consider the potential: 



v(q) = o + + x 3i3) , 



(9) 



where the Ai are the eigenvalues of the first-order deforma- 
tion tensor S^l = (f.ab', note the difference of sign from 
the definition in M95. They are ordered as follows: Ai > 
A2 > A3. Because of Poisson equation (^), 5i = S(to)/bo = 
Ai + A2 + A3 = const. Equation Q represents the potential 
of a homogeneous ellipsoid in its principal reference frame. 

With this potential it is possible to solve all the Poisson 
equations for the perturbing potentials. All the details of the 
calculations are contained Appendix B. The solutions can be 
directly found by using the so-called 'local forms' of the per- 
turbing potentials, which can be found in Buchert & Ehlers 
(1993), Buchert (1994) and Catelan (1995). Briefly, it is pos- 
sible to find vector fields, functions of the initial potential 
and its derivatives in t he point considered, which solve the 
spatial equations (A13), but are generally not irrotational. 
The local form S' 2 ^' of the second-order displacement is: 



where both qi and q2 are collapsed, relative to a scale L, at 
the time b > b c . This definition is very easy to implement in 
an N-body simulation. 



S w = [V^(VV)-(V<^-V)V^] + R W =S W +R W .(10) 

The vector R' 2 ', divergenceless, is added to the local form 
S' 21 '' in order to keep the vector S^ 2 ' irrotational; it contains 
the deep non-locality of gravitational dynamics. Analogous 
expressions for the 3a, 3b and 3c contributions are reported 
i n Ap pendix B. The local forms are solutions of equations 
( |A13| ) if they are irrotational by themselves; this happens 
only for a restricted class of initial conditions, to which the 
homogeneous ellipsoidal potential belongs. Then the local 
forms can be used to promptly obtain the perturbing poten- 
tials of the ellipsoid. 

The 2nd-order local contribution to the deformation 

(2L) 

tensor, S a b , contains terms with second derivatives of if 
and mixed terms with first and third derivatives. In the ho- 
mogeneous ellipsoidal case only the second derivative terms 
survive. The same happens for the 3aL and 3bL terms, while 
the 3cL term is null (see Appendix B). I call ellipsoidal parts 
those terms of the local contributions to the deformation 
tensor which survive in the ellipsoidal case: 
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M2E) 
D a,b 
c (3aE) 



^(36E) 
(3cE) 



c 



s. 



1 

2 
0; 



a(2E) , o(2E) 

S bc 'lfi,ac + <P,abS c , c 



(11) 



<p ab is the cofactor matrix of tp ia b- These terms can be con- 
sidered as a truncation of the local form, when all derivatives 
of ip greater than the second are neglected. (Note that 2nd- 
order local and ellipsoidal displacements, S 2L and S 2B , are 
equal, as the local form contains up to second derivatives of 
the initial potential; the differences appear in the deforma- 
tion tensor. Third-order local and ellipsoidal displacements 
are instead different.) 

If the very weak Q dependence of the time functions is 
neglected (see equations A.7), the collapse equation J = is 
simply an algebraic equation in the variable b(t) (the grow- 
ing mode), of order equal to the Lagrangian order. All the 
technical details are reported in Appendix B. In the spher- 
ical case, in which Ai = A2 = A3 = Si/3, the collapse times 
predicted by the Lagrangian series up to the 1st-, 2nd- and 
3rd-order, are b c Si = 3, 2.27, 2.05. The exact Q — 1 solution 
is b c — a c = 1.69; it can be appreciated that the results con- 
verge to the exact value: the difference is reduced from 1.31 
to .36. Other cosmologies lead to very similar values of b c , 
which differ from by not more than 3 per cent (Lilje 1992). 
For the Lagrangian perturbation scheme, spherical symme- 
try is the hardest to deal with, so we expect faster conver- 
gence in more general cases (see the discussion in M95). 

The first-order solution for general ellipsoids is simply 
60 = 1/Ai; it has been amply discussed in M95 (note that, 
due to the different signs of the Ai eigenvalues, b c was — I/A3 
in M95). The second-order solution for initially overdense 
ellipsoids is: 



;|; ,, _ 7Ai - ^/7Ai(Ai +6S1) 



u c 



(12) 



3Ai(Ai-<5i) 

If Si < 0, the second-order equation gives meaningful solu- 
tions only if Si > — A1/6, i.e. only for relatively small under- 
densities. Other solutions exist which make even a spherical 
void collapse! These had already been noted by Sahni & 
Shandarin (1996), and show how the 2nd-order Lagrangian 
scheme is unreliable for initially underdense perturbations. 
The reason for this behaviour is that 2nd-order perturba- 
tions, being of even order, do not properly recognize voids, 
making them collapse. 

The third-order solution is the smallest non-negative 
solution of the equation: 



A 1 6i 3) -^A 1 (^~A 1 )(6i 3) ) 2 



(13) 



V 126 + 84 



r \i6i{6i-\i))(b?>) 



«,(3h3 







where /13 = A1A2A3. Though it is not straightforward to 
choose analytically the right root of this equation, it is very 
easy to find it with a computer. The third-order equation 
gives meaningful solutions also for initially underdense el- 
ements, thus avoiding the problems encountered with the 
second-order solutions. This shows that third-order pertur- 



bation theory is necessary when predicting the instant of 
collapse of general mass elements. 

Both at second and third-order, the first axis to collapse 
is the one corresponding to the largest A eigenvalue, Ai ; this 
is by itself an indication of convergence, as it means that 
the first-order (Zel'dovich) approximation, which predicts 
collapse first on the 1-axis, makes the greatest contribution 
to collapse dynamics. 

The calculations of the exact collapse time of the el- 
lipsoid have been performed as in M95; the equations are 
reported in Appendix B. To improve the calculations, the in- 
tegration has been divided into two parts: after decoupling, 
defined as the instant when the density starts to grow, the in- 
tegration variable has been changed from (logarithm of) b to 
(logarithm of) the density; the integration has been stopped 
at 8 — exp(15). This had led to a modest but appreciable 
improvement in precision. The precision of the calculation 
depends on the initial conditions: the spherical collapse is 
predicted with a precision of around 0.2 per cent, while the 
'pancake' collapse, when A2 = A3 = 0, is recovered with 
a precision of 8 per cent. The calculations have been per- 
formed only for 8i = l or —1; the other cases can be found by 
using the scaling relation: 



MAi,A 2 ,A 3 ) = kb c (k\i,k\ 2 ,k\ 3 ) 



(14) 



Figs, la, b and c show the various b c curves of initially 
overdense and underdense ellipsoids. The x and y variables, 
which range from to 00, are defined as x = Ai— A2, y — A2 — 
A3; they give a measure of the initial shear of the ellipsoid 
(see M95). A number of things can be noted: 

(i) The predictions at increasing Lagrangian orders al- 
ways converge to the exact value (within the numerical er- 
rors quoted above) ; in the initial underdensity case, only the 
odd Lagrangian orders converge to the solutions. 

(ii) The convergence is very fast for large shears; in this 
case (for initially overdense ellipsoids) , the third-order solu- 
tion does not much improve the agreement with the numer- 
ical solution, with respect to the second-order one. 

(iii) Initially underdense ellipsoids can collapse if the 
shear is large enough; in this case the third-order predic- 
tion is always sufficiently accurate. 

These results can be used to predict the collapse time 
of a homogeneous ellipsoid in a fast and accurate way (prob- 
ably more accurate than numerical integration). It is neces- 
sary to correct the second- and third-order predictions in or- 
der to reproduce the correct behaviour in the quasi-spherical 
overdense cases. This correction can expressed in the follow- 
ing way: 

bi" c) = 



A exp(- 



by) 



(15) 



where n=2,3, and the three coefficients take on the following 
values: 



2nd ord 3rd ord 

A = 0.580 or 0.364 

a — 5.4 or 6.5 

b = 2.3 or 2.8 . 



(16) 



This corrections are applied only when Si > 0; no correction 
is applied when Si < 0. 
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Figure 1. Collapse times b c of homogeneous ellipsoids according to Lagrangian perturbation theory, compared with a numerical inte- 
gration; x = Ai — A2 and y = A2 — A3, (a): initial overdensity 5; = l.(b) and (c): initial underdensity 5; = — 1; 2nd-order predictions have 
been put in a separate figure for the sake of clarity, (d): open universe, initial overdensity 5; = 3; the Q dependence of the time functions 
is neglected. 



All these calculations have been repeated with two other 
background models, namely the open Qq = 0.2 model and 
the flat f2o = 0.2 model with cosmological constant. The La- 
grangian calculations are qui te insensitive to the cosmology, 
as shown by equations (A7) (note that the pancake case, 
x — Si and y — 0, being exactly predicted by Zel'dovich, is 
always independent of the cosmology). As expected, the nu- 
merical calculations give very similar results for the b c func- 
tion, with slightly larger errors when the background density 
becomes small and the ellipsoid takes a long time to collapse. 
Fig. Id shows the fi = 0.2 case; an initial overdensity Si = 3 
has been chosen to allow the collapse of all the ellipsoids 
reproduced in the figure. Then, the above-mentioned results 
for b c can be used in any of the cosmologies checked here. 

In conclusion, the Lagrangian series can be used to accu- 
rately approximate the collapse of a homogeneous ellipsoid. 
On the other hand, homogeneous ellipsoid collapse can be 
seen as a particular truncation of the Lagrangian series. Its 
ability to reproduce the collapse time of generic perturba- 



tions will be shown in the next section. In this case, ellip- 
soidal collapse has to be considered not as a description of 
the dynamics (or even the geometry!) of an extended region, 
but as an approximate description of the local dynamics of 
a mass element. 

Before going on, it is useful to comment on the concept 
of 'locality'. From a mathematical point of view, the evo- 
lution of a continuous system is local if its evolution equa- 
tions are ordinary differential equations, i.e. with no partial 
derivatives^]. Thus, any point's trajectory is fully determined 
by its initial conditions, i.e. it is not necessary to evolve 
all the trajectories together. From a physical point of view, 
non-locality can be associated to the long-range character 
of gravitational forces: the fate of a mass element depends 
on all the other mass elements. The non-local character of 
gravitational dynamics has recently been stressed again in a 



* It is assumed that there is no explicit coupling between field 
values at different points 
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paper by Kofman & Pogosyan (1995). Nonetheless, some dy- 
namical approximations predict 'local' dynamical evolution. 
The unrealistic spherical model is such a one: the fate of a 
spherical perturbation depends only on its initial density. At 
variance, the Zel'dovich approximation gives some informa- 
tion on non-local tidal forces: the motion of a mass element 
depends on all other elements. This information is contained 
in the gravitational potential, which is related to the den- 
sity by a non-local Poisson equation. Thus the Zel'dovich 
approximation is physically non-local. But, once the initial 
peculiar potential is known, the dynamical evolution is local, 
in the mathematical sense given above. In other words, the 
initial conditions contain non-local information, while the 
evolution is local. The same holds true for the homogeneous 
ellipsoidal collapse model. Analogous considerations can be 
applied to the Lagrangian perturbation theory at any or- 
der, as the equations for the perturbative terms are separa- 
ble in space and time, as recently demonstrated by Ehlers 
& Buchert (199 6). If this is the case, the non-local space 
equations (A15) give non-local initial conditions, while the 
subsequent dynamics is independently determined for every 
single trajectory. But, while it is relatively simple to ob- 
tain the statistical distribution of the initial conditions for 
the Zel'dovich approximation, namely of the Xi eigenvalues, 
the analytical determination of the statistical distribution 
of other contributions leads to discouraging mathematical 
difficulties. Nonetheless, as a consequence of the mathemat- 
ical 'locality' of Lagrangian perturbations, the perturbative 
terms have to be calculated just once at the beginning, mak- 
ing this kind of calculation dramatically faster than usual 
N-body simulations. 



5 COLLAPSE TIME IN THE GAUSSIAN 
FIELD CASE 

In this section, Gaussian fields with scale-free power spec- 
tra are considered. To proceed analytically, the equation 
J = det(5 a f, + S at b) = ought to be solved for general 
initial potentials, and the probability distribution function 
(hereafter PDF) of its smallest non-negative root ought to 
be obtained. As a matter of fact, it is very hard to get the 
whole PDF of the various contributions to the deformation 
tensor. Even neglecting the non diagonalizable 3c term, the 
various contributions to S a ,b are not diagonal in the same 
frame. It is definitely convenient to get the PDF of the col- 
lapse times by constructing realizations of Gaussian fields in 
cubic grids. 

In these calculations the grid does not need to be very 
large: what is needed is a sufficient degree of non-locality, 
which is provided even by small grids. 16 3 and 32 3 grids 
have been used, with identical results. Initial potentials have 
been simulated, following what is done for initial conditions 
of N-body simulations. Power spectra have been chosen as 
Pip(k) oc k n ~ 4 , with n=— 2, -1, and 1; these correspond to 
matter perturbation spectra Ps(k) oc k n . Spectra have been 
normalized with a = 1, where a is the total density variance. 
Of course, any other normalization can be obtained by a 
time rescaling. 

Ten realizations have bee n pe rformed for every power 
spectrum. Poisson equations (A13) have been solved with 
Fast Fourier Transform (FFT) techniques; FFT have also 



been used to calculate the derivatives of the various poten- 
tials. The growing mode b has been used as a time variable, 
and the Q dependence of the time functions b n (equations 
A7) has again been neglected. The following collapse time 
estimates have been calculated for every point: spherical 
collapse (hereafter SPH), lst-order or Zel'dovich approxi- 
mation (1ST), 2nd-order (2ND), 3rd-order (3RD) and el- 
lipsoidal collapse (ELL). SPH, 1ST and ELL collapse times 
have been calculated analytically, on the basis of the A eigen- 
values of (p t ab at any point; ELL has been calculated at third- 
order and corrected around the spherical value as in equation 
(E^). 2ND and 3RD have been calculated by looking for the 
instant at which J < 0, then using conventional root-finding 
algorithms. It is possible that J becomes negative and then 
positive soon after; these events can be lost if the search is 
not fine enough. As a matter of fact, a very small number of 
points, on the order of a few times 10~ 5 of the total number, 
were missed by the searching algorithm. 

Buchert, Melott & Weifi (1994) state that third-order 
terms are very sensitive to numerical calculations. The same 
thing has been found in these analyses; third-order calcula- 
tions are not expected to be very precise. This has been 
noted especially for the transversal 3c contribution. Any- 
way, that contribution has been found not to influence the 
collapse time appreciably, so it has been neglected in the 
calculations. 

From every set of 10 realizations, nearly 10000 points 
have been randomly extracted to analyze the statistics of 
collapse times. Figs. 2a-j and 3a-j show the scattergrams of 
the five collapse estimates, for n = —2 and n = 1. Points 
which are predicted not to collapse have been assigned a 
small negative collapse time, -0.1. Then, the points which 
are predicted to collapse according to one prediction, but 
not according to the other, are recognizable in the scatter- 
grams as horizontal or vertical rows of points. These points 
will be called discordant in the following. Mean values and 
dispersions around the bisector of non-discordant points are 
superimposed on the scattergram. Finally, Figs 2 and 3 fo- 
cus on the interesting zone b c < 3, which includes the points 
forming the large-mass part of the MF. 

Many conclusions can be drawn from Figs. 2 and 3: 

(i) As expected, SPH correlates with the other predic- 
tions for the fastest collapsing points (in agreement with 
Bernardeau 1994), but it badly overestimates the collapse 
time in general cases; moreover, many points (those with 
b~i < 0!) are incorrectly not predicted to collapse. Then, 
non-locality strongly accelerates collapse with respect to the 
spherical case, in line with the conclusions of M95. As a con- 
clusion, spherical collapse is not suitable, even statistically, 
for describing gravitational collapse. 

(ii) The 1ST - 2ND and the 2ND - 3RD correlations at 
small collapse times are increasingly good (though the for- 
mer has a considerable scatter); this demonstrates the con- 
vergence of the Lagrangian series in predicting the collapse 
time of the fast collapsing points. The 1ST - 3RD correlation 
is similar to that of 1ST - 2ND. 

(iii) The discordant points in the 2ND - 3RD scattergram 
are either some initially slightly non-negative ones, for which 
2ND does not find any solution, as in the ellipsoidal case, 
(and then be ND —— 0.1), or voids which are incorrectly pre- 
dicted to collapse by 2ND (and then b [ c 3RD) =~0.1). This 
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Figure 2. Scattergrams of the various collapse time estimates b c , for a Gaussian field with n = —2. 



shows the necessity of third order calculations for the col- 
lapse time. The same features are recognizable in the other 
2ND scattergrams. (N.B. The very small number of discor- 
dant points, less than 10 in 10000, which collapse according 
to 1ST and ELL, but not according to 3RD, are the ones 



missed by the algorithm which looks for J — 0; see the dis- 
cussion above.) 

(iv) 3RD accelerates the collapse of the points with b c > 1 
with respect to all the other predictions. Given the uncer- 
tainties connected with third-order calculations, this feature 
is not considered very robust. Moreover, 3RD predicts the 
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Figure 3. Scattergrams of the various collapse time estimates b c , for a Gaussian field with n = f . 



collapse of more points than 1ST. Again it is not clear at all if 
this feature, regarding points completely outside of the con- 
vergence range of the Lagrangian series, has some meaning. 
Probably Lagrangian perturbations are not a good means 
for determining the fraction of collapsed mass in a smooth 
universe. 



(v) ELL shows an encouraging correlation with 2ND and 
an even better one with 3RD, when b c is small. Note how 
1ST, 2ND and 3RD, when b c is small, converge to a solu- 
tion which tightly correlates with ELL. This has two implica- 
tions, namely that the Lagrangian series probably converges 
to the true solution and that ELL can be used as a realis- 
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tic estimate of the collapse time. Moreover, ELL tends to 
underestimate the collapse time for b c > 1, slightly with re- 
spect to 2ND and strongly with respect to 3RD. This would 
mean that the non-locality not contained in the ELL esti- 
mate speeds up the collapse. 

(vi) The results are essentially independent of n; the weak 
differences that can be visible in the scattergrams will be 
quantified in the following. 

These results are strictly correct for an Einstein de- 
Sitter Universe, in which case the scale factor can be used as 
a time variable, b c = a c . Due to the very weak ^-dependence 
of the time functions, when expressed in terms of b, the 
results for Q ^ 1 cosmologies are nearly indistinguishable. 
The most important effect, in the open case, is that the 
growing mode saturates at the value b = 5/2; a mass element 
with that collapse "time" collapses at an infinite physical 
time, and larger b c imply no collapse. 

As noted in Section 4, ELL can be considered as a trun- 
cation of the Lagrangian series, or, better, of its 'local forms'. 
It is interesting to check whether the use of the full local 
forms of the Lagrangian series improves the agreement of 
the approximation to the complete calculation. Here only 
the second-order terms will be considered; in this case, 2nd- 
order local and ellipsoidal displacements are equal, the dif- 
ferences come into the deformation tensor. The 2nd-order 
deformation tensor, its local and its ellipsoidal parts have 
been calculated for every point of a 16 3 realization with 
power spectrum n = 0. Fig. 4 shows the scattergrams of 
one of the diagonal elements of the 2ND deformation ten- 

(2) (2L) 

sor, Sf [ versus its local form, S{ 1 , and its ellipsoidal part, 

(2E) 

S\ 1 . Fig. 4 also shows the same pictures for one of the off- 
diagonal elements; other diagonal and off-diagonal elements 
obviously behave identically. The local diagonal elements 
correlate with the 2ND ones just slightly better better than 
the ellipsoidal ones, while no correlation of the off-diagonal 
terms is visible in any case. Fig. 5 shows the corresponding 
collapse times; here the ellipsoidal collapse time is given by 
equation (|l2j), without any correction. Again, the ellipsoidal 
part reproduces the full 2ND collapse time nearly as well as 
the local form. As a conclusion, it is convenient to truncate 
the local forms as in Section 4, as it leads to a great simpli- 
fication of the calculations and to a similar agreement with 
the full Lagrangian terms. 



5.1 The inverse collapse time PDF 

To quantify the PDF of the collapse times, it is convenient 
to consider their inverses, as they are better behaved: the 
inverse collapse times are large, but of order one, for fast 
collapsing points, and become smaller and smaller for slowly 
collapsing ones; the passage from collapse to non collapse is 
not at infinity, as for b c , but at 0. We define: 

F(q) = fe- 1 (q). (17) 

In the SPH case, F = 5;/1.69; in the 1ST case, F = Ai. 
This definition is also convenient for determining the MF 
(see paper II). 

Fig. 6 shows both the cumulative and differential F 
PDFs for n = —2 and 1; the cumulative curves give a 
binning-free picture of the PDFs, while the differential ones, 



in logarithmic scale, better exhibit the behaviour in the rare 
event tail. The following things can be noted: 

(i) The SPH curve is quite different from all the others, 
even at the high F tail: as in M95, a systematic departure 
from spherical collapse influences also the statistics of rare 
events. 

(ii) In the range F > 1, the 1ST, 2ND and 3RD curves 
show a monotonic shift toward large F values, which can 
be interpreted as convergence toward a solution. This is not 
true for F < 1; the bad behaviour of 2ND in initial under- 
densities is the probable cause. 

(iii) ELL and 3RD nearly coincide down to F = 1.2, 
which corresponds to b c — 0.83, and overall have a sim- 
ilar behaviour; ELL slightly makes more mass collapse at 
large F values, because 3RD slightly underestimates quasi- 
spherical collapses. The main differences come out in the 
range where the convergence of the Lagrangian series is not 
guaranteed, but both 1ST, 2ND and 3RD have larger medi- 
ans than ELL; so ELL probably underestimates the collapses 
around F = 0.5 or b c — 2. 

(iv) The n dependence can be appreciated in Fig. 6. As 
expected, SPH, 1ST and ELL are essentially independent 
of n (the \i distribution is fixed once the mass variance is 
fixed), while the n dependence of 2ND and 3RD is weak. 
Moreover, the difference between ELL and 3RD is smaller 
for smaller n, i.e. when more large-scale power is present, 
while we would expect the opposite if the difference of 3RD 
from ELL were due to non-locality induced from large scales. 
In the following the weak n-dependence of 3RD will be ne- 
glected. 

(v) The range F > 1, where the Lagrangian series con- 
verges to ELL, corresponds to more than 10 per cent of the 
points, i.e., of the mass. 20 per cent of the mass is found in 
F > 0.8, where ELL and 3RD are not very different, while 
all the medians can be found around F ~ 0.5. So the conver- 
gence of the Lagrangian series and its agreement with ELL 
take place for a significant quantity of mass, covering more 
than the large-mass tail of the MF, while the appreciable 
differences between ELL and 3RD affect the power-law part 
of the MF, which is plagued by a number of other problems 
(see paper II). 

To calculate the MF, an analytical expression for the 
F PDF is needed. To find it, it is useful to look for a co- 
ordinate transformation which maps the given PDF to a 
Gaussian one, with zero mean and unit variance. Consider 
then the quantity F with its PDF Pf(F), and a normal 
Gaussian Pg(x) = exp(— a; 2 /2)/v / 27r. If the function Pf is 
well-behaved enough, a transformation x(F) will exist such 
that: 

P F (F)dF = P G (x(F))dx(F) = -^ e - x(F)2/2 dx{F) . (18) 

V27T 

The function x(F) can be found as a solution of the ordinary 
differential equation: 

^)=v^P^)exp(^). (19) 

This equation can be easily solved with a computer. The 
initial condition of the equation is best set by integrating 
from the median of the Pf distribution, setting x(F) — at 
that point. 
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Figure 4. Scattergrams of diagonal (<p,ll) and non-diagonal (</>,i2) matrix elements of the 2nd-order contributions to the deformation 
tensor, calculated either exactly or by using local forms or their ellipsoidal parts. 



Equation ( J19| ) has been solved for the ELL and 3RD 
distributions. The 3RD distribution has been mediated over 
the four spectral indexes considered, so as to neglect the 
n dependence. The results are shown in Fig. 7a: in both 
cases the transformations for F > 1 are accurately linear, i.e. 
the large- F parts of the F distributions are Gaussians.(Note 
that the weak rise of the x(F) curve in the large- F end is 
an artifact). The x(F) transformations are accurately fit by 
the following expressions: 



x(F)ell = -0.69 + 1.82F - 0.4(erf(-7.5F + 1.75) + 1) (20) 
x{F)- ARD = -1.02 + 2.07F - 0.75(erf(-3F + 1.18) + 1) ; 



these analytical fits do not reproduce well the x(F) curves 
for F < 0.2; on the other hand, that part of the PDFs is 
very uncertain. The first two terms represent the linear fits, 
valid for large F values. Fig. 7b shows the two ELL and 
3RD PDFs, together with their Gaussian (coming from the 
linear transformations) and non-Gaussian (coming from the 
complete transformation) fits. Both non-Gaussian distribu- 
tions, compared with the Gaussian ones, show a peak at 
small F; the two peaks have similar heights but slightly dif- 
ferent positions. These correspond to the falling tail in the 
transformation curves x(F). 
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Figure 5. Scattergrams of 2ND collapse time estimates b c , calculated either exactly or by using local forms or their ellipsoidal parts. 



6 SUMMARY AND CONCLUSIONS 

In this paper, the dynamical part of a new theory of the mass 
function of cosmic structures has been presented. An accom- 
panying paper, paper II, develops the statistical machinery 
necessary to calculate the MF. The key point of all the dy- 
namical analysis is the use of a fiuidodynamical Lagrangian 
framework to predict the fate of a mass element of a smooth 
field. The Lagrangian perturbative scheme has been used 
to approximate the real dynamics of the smoothed cosmo- 
logical fluid up to orbit crossing. The strongest hypothesis 
of the whole theory is that small-scale structure does not 
influence the dynamics of larger collapsing scales. The OC 
instant has been chosen as suitable definition of collapse of 
a mass element. This definition has been amply discussed. 

To check the performances of Lagrangian perturbations 
in predicting the OC instant, when all the terms of the se- 
ries are of the same order, the collapse of a homogeneous 
ellipsoid has been considered. The Lagrangian series con- 
verges fast in predicting the collapse of a homogeneous el- 
lipsoid, except for quasi-spherical ones; a simple correction 
can be applied in this case. Third-order terms are necessary 
to give good collapse estimates of initial underdensities. On 
the other hand, ellipsoidal collapse can be considered as a 
truncation of the whole Lagrangian series. In this case, as in 
M95, ellipsoidal collapse has to be considered an approxima- 
tion of the local collapse of a point mass, not of an extended 
region of approximate ellipsoidal shape. 

The collapse of scale-free Gaussian fields simulated on 
a 32 3 grid has been considered. The Lagrangian series has 
been found to clearly converge to a solution for the fast- 
collapsing points, about 10 per cent of the total mass. Ap- 
proximate convergence has been observed for about 50 per 
cent of the mass. The homogeneous ellipsoidal model has 
been found to strongly correlate with the third-order La- 
grangian collapse in the same range of validity. The ellip- 
soidal model can thus be used as a fast and easy to im- 
plement approximation of collapse dynamics. The spherical 



model, instead, has been found to behave differently from 
the other predictions. The inverse collapse-time distribu- 
tions have been calculated for the third-order Lagrangian 
and ellipsoidal predictions. These will be necessary for de- 
termining the MF. 

With this dynamical description of collapse, it is pos- 
sible to construct a MF, fully based on realistic dynamics. 
The role of the initial density is taken by the inverse collapse 
time F. This quantity is by no means a Gaussian process, 
but is a complicated non-linear functional of a parent Gaus- 
sian process. This introduces a number of complications in 
the statistical treatment of the MF, which will be faced in 
paper II. 

The following points, raised by the analysis presented 
here, are relevant when constructing a MF: 

(i) As the main quantity, F, is (the inverse of) a time, 
any threshold F c in this theory simply specifies the time at 
which the MF is examined; there is no free S c parameter. 

(ii) The dynamical predictions are strictly punctual; in 
other words, a point collapses if it is predicted to collapse 
(at a given scale), not if a neighboring collapsing point is 
able to involve it in its collapse. 

(iii) Smoothing is necessary because of the truncated na- 
ture of the dynamical approximations used. Thus the shape 
of the filter has to be chosen in order to optimize the per- 
formances of the dynamical predictions; usually Gaussian 
filters are suggested. 

The collapse time, which is needed to determine the 
MF, is not the only information one can obtain on a col- 
lapsing mass element. In practice, all the trajectories of the 
collapsing points are known. This has the consequence that 
much more than a MF can be calculated. A possible ap- 
plication is the determination of the angular momentum of 
a collapsing region; Lagrangian perturbation theory is an 
appropriate framework for estimating such a quantity, as 
Catelan & Theuns (1996a,b) have recently shown. The rich- 
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Figure 6. Cumulative and differential PDF of various inverse collapse time estimates F, for n = —2 or n = 1. 



ness of the dynamical information makes it possible to put 
further dynamical constraints on the collapsing regions, if 
special classes of objects are required. 

To finally decide on the validity of the MF dynamics 
presented here, the predictions of this theory have to be 
tested against N-body simulations. As collapsed regions are 
defined to coincide with orbit-crossed regions, they have to 
be sought for by a suitable algorithm; equation (|^) could 
provide such an algorithm. This MF is expected to carefully 
reproduce the N-body MF as long as Lagrangian perturba- 
tion theory is expected to give a correct description of col- 
lapsing dynamics, i.e. for large masses, comprising at least 
10-20 per cent of mass, and when few small-scale structures 
are present, i.e. with small spectral indexes. 
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APPENDIX A: LAGRANGIAN 
PERTURBATIONS 

This appendix contains a number of technical points on La- 
grangian perturbation theory. 



Schechter P., 1974, ApJ, 187, 425 (PS) 
, Teukolsky S.A., 1990, Computers in Physics Jan/- 
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Peculiar velocity, acceleration and density contrast of a 
mass element can be written in terms of the S field as: 



v(q,i) = a(t)dS(q,t)/dt 
g(q,t) = a{t)d 2 S(q,t)/dt 2 +2tfv(q,t) 
l + *(q,t) = J(q,t) _1 [l + *(q.*o)]i 



(Al) 



here v(q, t) is the peculiar velocity of the element q, g(q,t) 
is its peculiar acceleration, 5 its density, d/dt denotes to- 
tal (Lagrangian) time derivative, a(t) is the scale factor of 
the background cosmology, H — a~ 1 da/dt is the Hubble 
parameter, to is an initial time and J(q, t) is the Jacobian 
determinant, equation (^). 

The evolution equations for the displacement field S can 
be written, following Catelan (1995), as: 

[(1 + V ■ S)S bd - S M + S£ d ]S a , d = a(r)[J - 1] (A2) 
Ea6o[(l + V ■ S)S bd - S b , d + S£ d ]S c , d = , (A3) 

where the dot denotes the Lagrangian derivative with re- 
spect to the time variable r = t~ 1 ^ 3 = a -2 , if flo = 1, or 

-1/2 



t = li-ai 



(A4) 



otherwise. e a tc is the Levi-Civita antisymmetric tensor, 
V = d/dq, S c is the cofactor matrix of S and the func- 
tion q(t) = 6/(t 2 + k) (k = — 1, or 1 for open, flat 



and closed models). The first equation (4.2) is an evolu- 
tion equation for S, i.e. the equati on o f motion of a mass 
element, while the second equation (A3) is the irrotational- 
ity condition for the peculiar velocity v in Eulerian space. 
The irrotationality condition restricts the set of solutions to 
the irrotational ones, which is reasonable in our cosmolog- 
ical context as any rotational mode is severely damped in 
the early linear evolution; see Buchert (1992) for a detailed 
discussion of this point. Another useful restriction on the so- 
lutions of equations (AS) and (A.3) is the initial parallelism 
of peculiar velocity and acceleration, which is also supported 
by the growing modes in the linear and quasi-linear regime 
(Buchert 1992; Buchert & Ehlers 1993). 

The perturbative expansion is performed as follows. We 
can write the displacement field S as: 

S(q,t) = eS (1) (q,t)+e 2 S (2) (q,t)+e 3 S (3) (q,t)+e)(e 4 ) ,(A5) 

where e is a small parameter. Putting this expression into 
the evolution equations ( AS , Ag] ) for S, and considering terms 



of order e, e and e separately, one finds equations for the 
various S'™' terms. It turns out that, at any order (as re- 
cently demonstrated by Ehlers & Buchert 1996), the solu- 
tions are separable in time and space. At first and second 
orders, the solutions are irrotational in Lagrangian space, 



i.e. the matrices and S K ^' b are symmetric. At the third 
order, the equation for S^ 3 ' is not separable as it stands, 
but it is possible to divide the S 1 - 3 ' term into three different 
modes, all separable in space and time: 

S (3) (q, t) = S (3a) (q, t) + S (3b) (q, t) + S (3c) (q, t) . (A6) 

The third term, not reported by Bouchet et al. (1995) who 
consider only those terms which contribute to S a ,a, is a 
purely rotational mode in Lagrangian space (S a 3 £ is an- 
tisymmetric), and its existence is necessary to guarantee ir- 
rotationality in Eulerian space. This fact can be understood 



? (2) 



in this way: the Lagrangian to Eulerian transformation is 
in general a non-Galileian one, so the rotational mode can 
be seen as an effect of inertial forces (see the discussions in 
Buchert 1994 and Catelan 1995). 

Perturbative terms are listed in the following. The equa- 
tions for the time functions b n contain both growing and 
decaying modes, which have to be consistently considered 
in the calculations; however, for present purposes only the 
growing modes for each b n are needed. The time functions 
are accurately described by the following expressions (exact 
for the first order): 

bi = -b(t) 

b 2 = -^bltt- a 

b 3 b = ~"^ b ^ C 
&3c = jrblCl~ d , 

where b(t) is the linear growing mode. In an Einstein-de 
Sitter background it is simply: 



(A7) 



b(t) = a(t) 



(A8) 



In an open Universe, it is convenient to express the growing 
mode in terms of the time variable r, equation ( A4): 

5 /, . „, 2 ,n /„ . t , (r — V 
2 



6(r) 



(l + 3(r 2 



1 + - In , 

2 Vr + 1 



(A9) 



The growing mode has been normalized so as to give b(t) ~ 
a{t) at early times. In a flat Universe with cosmological con- 
stant, it is useful to use h — coth(3iifoi\/l — Ho/2) as time 
variable: 



b(h) =h (x 2 {x 2 ~ l) 173 ) -1 ^ 



(A10) 



The coefficients a, b and c in equation (A.7) have been 
calculated in Bouchet et al. (1992) and Bouchet et al. (1995), 
and are a=2/63, b—4/77 and c=2/35 in the non-flat cases, 
0=1/143, fe=4/275 and c=269/17875 in the flat A / 
cases. The above-cited authors have not estimated the in- 
dependence of the 3c term, as they do not take that term 
into account. However, as clarified in the text, the 3c term 
can be safely neglected. It can be appreciated that the b^lb 2 
and b'i/b 3 terms weakly depend on Q when Q ~ 1; then, in a 
Universe with fio > 0.2, as our Universe appears to be, the 
Q dependence in equations ( [A7| ) can be safely neglected. 

The spatial equations for the S' n ^(q) terms are Poisson 
equations. It is convenient to express S* 1 ', S*- 2 ', S' 3 "- 1 and 
S' 36 ' in terms of scalar potentials, and s' 311 ' in terms of a 
vector potential: 



s (n) = y^tn)^ n=1]2 ,3a,36 
S (3c) = V x v> (3c) . 



(All) 



Defining the principal and mixed invariants of one or two 
tensors as follows: 



m(Aab) = tr(A a b) = A aa 
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fJ-2{A ab , B ab ) = -(A aa B bb - A ab B ab ) (A12) 

H2{A ab ) = fJ-2(A ab , A ab ) 
^3{A ab ) = det(A ab ) 

(note that fj,i((p, ab ) = V 2 <^), it is possible to write the fol- 
lowing equations for the potentials: 



(2) 



vv 3a) 
vV 36) 

vVi 3c) 



= <p 

(1) (2) 
= Safeco c > d 6 ■ 



(A13) 



The first equality is a consequence of initial conditions. Cate- 
lan (1995) gives expressions for the Fourier transforms of the 
solutions of all these equations, which are useful for calcu- 
lating, e.g., mean values of the perturbing terms. 



APPENDIX B: ELLIPSOIDAL COLLAPSE 

Let y(q) be a quadratic potential in its principal reference 
frame: 



y(q) = o ( Xiq i + X ' 2q 2 + *3?f) 



(Bl) 



It is easy to calculate all the perturbative terms in this 
case, especially if the local forms are used (see Buchert & 
Ehlers 1993, Buchert 1994 and Catelan 1995 for further de- 



tails) : 








(2L) 
<P,a 


— (P,a<P, bb — <P, 


a b y>, b 






c 

= f.ab'fi.b 




(B2) 






(2) , (2) 
- V,bf ,ab + <P)a <P,bb 


(2) \ 
V,b <P,ab) 




= V2fe^S 


(2) \ 





<p ab is the cofactor matrix of <p : ab- These local parts are exact 
solutions in our case, as they are irrotational; their expres- 
sions can be considerably simplified, as all the derivatives 
beyond the second vanish. The outcoming contributions to 
the deformation tensor have been given in equation (|ll|). 
A technical remark on the 3b contribution: Buchert (1994) 
divides this contribution into two parts, weighted by two co- 
efficients whose sum is equal to one (see his equations 27). 
These two coefficients could be varied in order to make the 
whole contribution irrotational. On the other hand, Cate- 
lan gives an expression analogous to the one given here (see 
his equation 44); this would correspond to a choice of 1/2 
for the two coefficients of Buchert. As a matter of fact, the 
two parts identified by Buchert are irrotational by them- 
selves in the ellipsoidal case, and have the same divergence, 
but are nonetheless different. Using only one or another, 
which would correspond to setting one coefficient to 1 and 
the other to 0, makes the Lagrangian series no to converge 
any more to the numerical solution. So, the choice of 1/2 for 
both coefficients seems the right one for ellipsoidal collapse. 

All the contributions to the deformation tensor are di- 
agonal in the same frame; the 3c contribution obviously van- 
ishes. The diagonal (1,1) components are: 



<p,u = 

(2) 

¥>,ii = 

(3b) 

¥>n = 



Ai 

Ai(A 2 + A 3 ) 
A1A2A3 

A1A2A3 + Ai<5; (A2 + A 3 )/2 



(B3) 



where 5; = Ai + A2 + A3. The Jacobian determinant van- 
ishes when one of its eigenvalues vanishes. Then, if the SI 
dependence of the time functions is neglected, it is possi- 
ble to write the J = equation as a third-order algebraic 
equation: 



1 - A,6 C - jgXi(5i - Xi)bt - + ^MSi - A*)) bl 



: 



(B4) 



where b is the linear growing mode, and ^3 = A1A2A3. Note 
how higher-order coefficients become increasingly smaller. 
In the spherical case, the equation reduces to: 



1 1 93 



(B5) 



which, if truncated at first, second, or third-order, gives the 
solutions: 



b^ = 3A 
b[ 2) = 2.27/5, 
b { c 3) = 2.05/5; 



(B6) 



The complete equation at first-order gives the well 
known Zel'dovich approximation: 



b c = 1/Xi 



(B7) 



It is apparent that the 1-axis, corresponding to the largest A 
eigenvalue, is the first to collapse. The second-order solution 
is: 



6 (2) 



TAi - V7Ai(Ai +65Q 
3Ai(Ai -5;) 



(B8) 



This solution is limited to 5; > — A1/6, i.e. to overdensities 
or weak underdensities. The other solution of the equation 
(with a plus in front of the square root) is either negative 
or larger than the chosen one, except when 5; < A^ < 0. In 
this case, which includes spherical voids, the second solution 
incorrectly predicts collapse. The bad behaviour of second- 
order perturbations in predicting the dynamics of voids had 
already been noted by Sahni & Shandarin (1995). Finally, 
it is possible to verify, by differentiating equation (B8) with 
respect to the Ai parameter, that the 1-axis is the first to 
collapse. 

To obtain the third-order solution, let y — 5i/Xi, and 
D = a*3A 3 /126 + 5y{y - l)/84; then let 

Q = (3(y - l) 2 - 196L»)/588_D 2 



R 



(2(y - l) 3 - 196(2/ - - 2744L» 2 )/5488L» 3 



When R 2 — Q 1 > 0, which is valid for spherical and quasi 
spherical perturbations, the solution is: 
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Q U/B? - Q 3 + \R\ 



1/3 

3 ^m\) + 



(B9) 



(y-i) 

14£> 



In this case, it is possible to show analytically that the 1- 
axis is the first to collapse. Otherwise, the solution has to 
be chosen as the smallest non-negative one between: 



(BIO) 



b c i = -2^/<2cos(6>/3) - (y- 1)/14D 
b c2 = -2 V / Qcos((6l + 27r)/3)- (y - l)/UD 
bcs = -2 v /Qcos((6» + 47r/3)-(2/-l)/14D ] 

where 9 = arccos(i?/ Ak /Q 3 ). It has been checked with the 
computer that the 1-axis is the first to collapse. 

The evolution equations of the homogeneous ellipsoid, 
which have been numerically integrated, are analogous to 
the ones in M95: let Vi(t) = a,i(t)qi be the physical i-th co- 
ordinate of the outer surface of the ellipsoid, in the principal- 
axes system; qt is the initial - Lagrangian - position and a^ 
is the expansion factor for the i-th axis. The universal scale 
factor is called a, without pedix. Then the evolution equa- 
tions for the a,i factors are: 



(2a 2 (l + (fi c 7 1 - l)a)) -1 ai 
in the open case and 



(Bll) 



1 5 b'i . , , 

— I h — 6 + X vi 

3 3 2 



d 2 a,i _ 1 - 2(^Q 1 - l)a 3 dm 
da 2 2o(l + (Qq 1 - l)o) da 

(2a 2 (l + (n f 7 1 -l)a))- 1 a l 



(B12) 



1 8 b[ ; , ,/ 

— I h — 5 + X vi 

3 3 2 



in the flat case with cosmological constant. Here (ao is the 
initial scale factor): 

5 = — 1 (B13) 

aia 2 a 3 

2 

b'i = -[aiajakRD(ai,a?,al) - 1] ij^j^k (B14) 



Kji — f 77 — doAi 

ao \3 



R D (x,y,z) = - 



dr 



2 Jo (r + x) 1 /2( T + y )i/2( r + 2 )3/2 



(B15) 
(B16) 



This Carlson's elliptical integral has been calculated by 
means of the routine given by Press & Teukolsky (1990). 
The initial conditions are: 



ai(ao) = a (l — aoK) 

-j^(ao) = — (oi(ao) — ao^O 
aa ao 



(B17) 
(B18) 



This paper has been produced using the Royal Astronomical 
Society/Blackwell Science I^TfrjX style file. 
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